
library(tidyverse)
library(gridExtra)

ABS_DIFF = function(df, ctr, i) {
  
  diff = df
  diff[,i] = df[,i] - ctr[,i]
  return(diff)
  
}

REL_DIFF = function(df, ctr, i) {
  
  diff = df
  diff[,i] = (df[,i] - ctr[,i]) / ctr[,i] * 100
  return(diff)
  
}

thin_wy = read_csv("../calibration/thinwyres.csv")
thin_bwy = read_csv("../calibration/thinwybres.csv")
thin_zwy = read_csv("../calibration/thinwyzres.csv")

thinN_wy = read_csv("../calibration/thinwyresN.csv")
thinN_bwy = read_csv("../calibration/thinwybresN.csv")
thinN_zwy = read_csv("../calibration/thinwyzresN.csv")

thinH_wy = read_csv("../calibration/thinwyresH.csv")
thinH_bwy = read_csv("../calibration/thinwybresH.csv")
thinH_zwy = read_csv("../calibration/thinwyzresH.csv")

thinHN_wy = read_csv("../calibration/thinwyresHN.csv")
thinHN_bwy = read_csv("../calibration/thinwybresHN.csv")
thinHN_zwy = read_csv("../calibration/thinwyzresHN.csv")

thin_wy = rbind(thin_wy, thinN_wy, thinH_wy, thinHN_wy)
thin_bwy = rbind(thin_bwy, thinN_bwy, thinH_bwy, thinHN_bwy)
thin_zwy = rbind(thin_zwy, thinN_zwy, thinH_zwy, thinHN_zwy)

#rm(thinN_wy, thinN_bwy, thinN_zwy, thinH_wy, thinH_bwy, thinH_zwy, thinHN_wy, thinHN_bwy, thinHN_zwy)


thin_bwy$T_sum = thin_bwy$T_sum / 1000
thin_bwy$E_sum = thin_bwy$E_sum / 1000
thin_bwy$ET_sum = thin_bwy$ET_sum / 1000
thin_bwy$dSD = -1 * thin_bwy$dSD

# streamflow metrics
thin_bwy$thin = if_else(thin_bwy$loc!="Control", paste(thin_bwy$loc, thin_bwy$int, sep=" "), thin_bwy$loc)
thin_bwy$loc = factor(thin_bwy$loc, levels=c("Control", "Thin Upslope", "Thin Riparian", "Thin Both", "Burn Upslope", "Burn Riparian", "Burn Both"))
thin_bwy$elev = factor(thin_bwy$elev, levels=c("Lo","Hi"), labels=c("Low-Elevation","High-Elevation"))


# calculate differences to control
compsf_list = split(thin_bwy, f=thin_bwy$thin)

compsf_diff = lapply(compsf_list, REL_DIFF, compsf_list[["Control"]], c(2:5,7:9,11,12))
compsf_absdiff = lapply(compsf_list, ABS_DIFF, compsf_list[["Control"]], c(2:5,7:9,11,12))
metsf_diff = do.call(rbind, compsf_diff) %>% filter(thin!="Control")
metsf_absdiff = do.call(rbind, compsf_absdiff) %>% filter(thin!="Control")

# streamflow time-series

ggplot(filter(metsf_diff, elev=="Low-Elevation", type=="Thin"))+geom_boxplot(aes(as.factor(t_wy), Q_sum, group=interaction(t_wy, int), fill=loc, alpha=as.factor(int)), outlier.size=.3)+
  geom_hline(yintercept=0)+facet_wrap(~loc, ncol=1)+scale_fill_brewer(palette="Dark2")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y="Annual Flow (% Deviation)", alpha="Intensity")+guides(fill="none")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot(filter(metsf_absdiff, elev=="High-Elevation", type=="Thin"))+geom_boxplot(aes(as.factor(t_wy), Q_sum, group=interaction(t_wy, int), fill=loc, alpha=as.factor(int)), outlier.size=.3)+
  geom_hline(yintercept=0)+facet_wrap(~loc, ncol=1)+
  theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Year since Treatment", y="Annual Flow (% Deviation)", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


ggplot(filter(metsf_diff, elev=="Low-Elevation", type=="Thin"))+geom_boxplot(aes(as.factor(t_wy), Q_max, group=interaction(t_wy, int), fill=loc, alpha=as.factor(int)), outlier.size=.3)+
  geom_hline(yintercept=0)+facet_wrap(~loc, ncol=1)+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y="Peak Flow (% Deviation)", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot(filter(metsf_diff, elev=="Low-Elevation", type=="Thin"))+geom_boxplot(aes(as.factor(t_wy), Q_min, group=interaction(t_wy, int), fill=loc, alpha=as.factor(int)), outlier.size=.3)+
  geom_hline(yintercept=0)+facet_wrap(~loc, ncol=1)+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y="7-day Minimum Flow (% Deviation)", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))


ggplot(metsf_absdiff)+geom_boxplot(aes(as.factor(t_wy), Q_sum, group=interaction(t_wy, int), fill=loc, alpha=as.factor(int)))+
  geom_hline(yintercept=0)+facet_wrap(~loc, ncol=1)+
  theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Years since Thinning", y="Annual Flow (Abs. Deviation)", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


ggplot(filter(metsf_diff, t_wy>2), aes(aspect, Q_min, group=interaction(aspect, int), fill=loc, alpha=as.factor(int)))+geom_boxplot()+facet_wrap(~loc)+
  theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Aspect", y="7-day Minimum Flow (% Deviation)", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


ggplot(filter(metsf_diff, t_wy>2), aes(elev, Q_min, group=interaction(elev, aspect), fill=loc, alpha=aspect))+geom_boxplot(outlier.size=.3)+facet_wrap(~loc)+
  theme_bw()+theme(legend.position=c(.9,.15))+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Elevation", y="7-day Minimum Flow (% Deviation)", alpha="Aspect")+ylim(-10,11)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


h = ggplot(filter(metsf_diff, elev=="High-Elevation"), aes(Q_max, Q_min, color=SWE_pk))+geom_point(size=1)+facet_wrap(~elev)+theme_bw()+theme(legend.position="bottom")+labs(x="% Change in Peak Flow", y="", color="Peak SWE (m)")
l = ggplot(filter(metsf_diff, elev=="Low-Elevation"), aes(Q_max, Q_min, color=SWE_pk))+geom_point(size=1)+facet_wrap(~elev)+theme_bw()+theme(legend.position="bottom")+labs(x="% Change in Peak Flow", y="% Change in Minimum Flow", color="Peak SWE (m)")
grid.arrange(l,h,ncol=2)

# compare water balance

metsf_wy = metsf_absdiff %>% select(-c(Q_min, Q_max, Q_CM, P_sum, SWE_pk, SWE_pre, SWE_pre2, melt)) %>% 
  group_by(t_wy, elev, int, loc, type, thin) %>% summarize_all(list(mean)) %>% 
  select(-c(scen, param, clim, aspect))

metsf_wb = gather(metsf_wy, key="Flux", value="Value", Q_sum, T_sum, E_sum, dSD)

ggplot(filter(metsf_wb, elev=="Low-Elevation", int==60, type=="Thin"), aes(t_wy, Value, fill=Flux))+
  geom_area()+facet_wrap(~loc, ncol=1)

ggplot(filter(metsf_wb, elev=="Low-Elevation", int==60, type=="Thin"), aes(t_wy, Value*1000, fill=Flux))+
  geom_col()+facet_wrap(~loc, scales="free_y", ncol=1)+theme_bw()+
  theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y=expression(Delta~"Water Balance Flux (mm/year)"))+
  scale_fill_manual(values=c("#E41A1C", "#984EA3", "#377EB8", "#4DAF4A"),
                    labels=c("Soil Storage","Evaporation","Streamflow","Transpiration"))


# tradeoffs

ggplot(data=filter(metsf_absdiff, type=="Thin"))+geom_hline(yintercept=0)+geom_vline(xintercept=0)+
  geom_point(aes(ET_sum*1000, Q_sum*1000, color=dSD*1000), size=1)+facet_grid(elev~loc)+
  geom_abline(slope=-1, intercept=0, linetype="dashed")+scale_color_viridis_c()+
  labs(x=expression(Delta~"ET (mm/year)"), y=expression(Delta~"Q (mm/year)"), color=expression(Delta~"S (mm/year)"))+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), 
                   strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))

# climate comparisons

ggplot(filter(metsf_diff, t_wy<3, elev=="Lo"), aes(P_sum, Q_sum, color=T_sum))+geom_point()+facet_grid(int~loc)+
  geom_hline(yintercept=0)+labs(x="Annual P (m)", y="Annual Q (% Deviation)", color="Change\nin Trans")

ggplot(filter(metsf_absdiff, t_wy<10, elev=="Low-Elevation", type=="Thin"), aes(P_sum, Q_sum*1000, color=SWE_pk))+geom_point()+facet_grid(int~loc)+
  theme_bw()+theme(axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  geom_hline(yintercept=0)+labs(x="Annual P (m)", y=expression(Delta~"Q (mm)"), color="SWE (m)")

ggplot(filter(metsf_absdiff, t_wy<3, elev=="Lo"), aes(melt, Q_sum*1000, color=T_sum))+geom_point()+facet_grid(int~loc)+
  geom_hline(yintercept=0)+labs(x="Snowmelt Day", y="Annual Q (Abs. Deviation, mm)", color="Change in\nTrans (mm)")

l=ggplot(filter(metsf_absdiff, elev=="Low-Elevation", type=="Thin"), aes(P_sum, Q_min, color=SWE_pre2))+geom_point()+facet_grid(elev~loc, scales="free")+
  theme_bw()+theme(axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  geom_hline(yintercept=0)+labs(x="Annual P (m)", y=expression(Delta~"Q min (mm)"), color="SWE (m)")

h=ggplot(filter(metsf_absdiff, elev=="High-Elevation", type=="Thin"), aes(P_sum, Q_min, color=SWE_pre2))+geom_point()+facet_grid(elev~loc, scales="free")+
  theme_bw()+theme(axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  geom_hline(yintercept=0)+labs(x="Annual P (m)", y=expression(Delta~"Q min (mm)"), color="SWE (m)")

grid.arrange(l,h)

ggplot(filter(metsf_absdiff, t_wy==1, type=="Thin"), aes(P_sum, Q_max, color=SWE_pk))+geom_point()+facet_grid(elev~loc)+
  theme_bw()+theme(axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  geom_hline(yintercept=0)+labs(x="Annual P (m)", y=expression(Delta~"Q (m)"), color="SWE (m)")


# parameter uncertainty
thin_bwy_p = thin_bwy %>% select(-thin, -type) %>% group_by(param, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% filter(elev=="Low-Elevation")
CVp = apply(thin_bwy_p[,7:9], 2, function(x){sd(x)/mean(x)})
CVp = data.frame(Var=colnames(thin_bwy_p[,7:9]), CVp)

metsf_diff_p = metsf_diff %>% select(-thin, -type) %>% group_by(param, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% filter(elev=="Low-Elevation")
CVp2 = apply(metsf_diff_p[,7:9], 2, function(x){sd(x)/mean(x)})
CVp2 = data.frame(Var=colnames(metsf_diff_p[,7:9]), CVp2)

# climate uncertainty
thin_bwy_c = thin_bwy %>% select(-thin, -type) %>% group_by(clim, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% filter(elev=="Low-Elevation")
CVc = apply(thin_bwy_c[,7:9], 2, function(x){sd(x)/mean(x)})
CVb = cbind(CVp, data.frame(CVc))

metsf_diff_c = metsf_diff %>% select(-thin, -type) %>% group_by(clim, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% filter(elev=="Low-Elevation")
CVc2 = apply(metsf_diff_p[,7:9], 2, function(x){sd(x)/mean(x)})
CVb2 = cbind(CVp2, data.frame(CVc2))


CVb_comp = rbind(mutate(thin_bwy_p, comp="Param"), mutate(thin_bwy_c, comp="Clim"))
CVb_long = gather(CVb_comp, key="Q", value="Flow", Q_max, Q_min, Q_sum, Q_CM) %>% select(-P_sum, -T_sum)
ggplot(CVb_long, aes(comp, Flow, fill=comp))+geom_boxplot()+facet_wrap(~Q, scales="free_y")+
  labs(x="Uncertainty", fill="Uncertainty")+scale_fill_brewer(palette="Accent")

CVb_comp2 = rbind(mutate(metsf_diff_p, comp="Param"), mutate(metsf_diff_c, comp="Clim"))
CVb_long2 = gather(CVb_comp2, key="Q", value="Flow", Q_max, Q_min, Q_sum) %>% select(-P_sum, -T_sum)
CVb_long2$Q2 = with(CVb_long2, ifelse(Q=="Q_max", "Peak Q", ifelse(Q=="Q_sum", "Annual Q", "Min Q")))
ggplot(CVb_long2, aes(comp, Flow, fill=comp))+geom_boxplot(outlier.size=.3)+facet_wrap(~Q2, scales="free_y")+
  labs(x="Uncertainty", y=expression("%"*Delta~Q))+scale_fill_brewer(palette="Accent", direction=-1)+theme_bw()+
  theme(legend.position="none", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  ggtitle("Uncertainty in Streamflow Response to Thinning")

# tables

metsf_diff$loc8 = with(metsf_diff, ifelse(grepl("Both",loc),"Both", ifelse(grepl("Upslope",loc),"Upslope", "Riparian")))
metsf_diff2 = metsf_diff %>% filter(elev=="Low-Elevation") %>% group_by(int, loc8) %>% 
  summarize_at(vars(Q_sum, Q_min, Q_max, Q_CM), list(mean=mean, sd=sd))


# stratum metrics

thin_wy$strata = thin_wy$stratumID
thin_wy$strata = factor(thin_wy$strata, levels=c(36815011, 36815021, 36815031, 36815041, 36815051, 38571011, 38571021, 38571031, 38571041, 38571051,
                                                 46901011, 46901021, 46901031, 46901041, 46901051, 48536011, 48536021, 48536031, 48536041, 48536051),
                        labels=c("Upslope Old-Growth", "Upslope Medium", "Upslope Thinned", "Upslope Shrub", "Upslope Grass",
                                         "Riparian Old-Growth", "Riparian Medium", "Riparian Thinned", "Riparian Shrub", "Riparian Grass",
                                         "Upslope Old-Growth","Upslope Medium", "Upslope Thinned", "Upslope Shrub", "Upslope Grass",
                                         "Riparian Old-Growth", "Riparian Medium", "Riparian Thinned", "Riparian Shrub", "Riparian Grass"))
thin_wy$thin = if_else(thin_wy$loc!="Control", paste(thin_wy$loc, thin_wy$int, sep=" "), thin_wy$loc)
thin_wy$loc = factor(thin_wy$loc, levels=c("Control", "Thin Upslope", "Thin Riparian", "Thin Both", "Burn Upslope", "Burn Riparian", "Burn Both"))
thin_wy$Site = factor(thin_wy$Site, levels=c("U","R"), labels=c("Upslope","Riparian"))
thin_wy$elev = factor(thin_wy$elev, levels=c("Lo","Hi"), labels=c("Low-Elevation","High-Elevation"))
thin_wy$patch = substr(thin_wy$stratumID, 6, 9)
thin_wy$sddiff = -1 * thin_wy$sddiff


# forest stand scale metrics
thin_wy$area = with(thin_wy, case_when(substr(stratumID,6,9)=="011"~.1, substr(stratumID,6,9)=="021"~.3, substr(stratumID,6,9)=="031"~.3, substr(stratumID,6,9)=="041"~.15, substr(stratumID,6,9)=="051"~.15))
thin_swy = thin_wy %>% group_by(Site, t_wy, param, clim, aspect, elev, int, loc, thin) %>% 
  summarize(standc = plantc %*% area, lai = lai %*% area, trans = trans %*% area, evap = evap %*% area, ET = ET %*% area, npp = npp %*% area, swe = swe %*% area, pcp = pcp %*% area, sddiff = sddiff %*% area, recharge = recharge %*% area, .groups="keep")
#thin_swy$thin = if_else(thin_swy$loc!="Control", paste(thin_swy$loc, thin_swy$int, sep=" "), thin_swy$loc)

ggplot(filter(thin_swy, loc %in% c("Control","Thin Both")), aes(as.factor(t_wy), lai, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.5)+geom_hline(yintercept=0)+facet_wrap(elev~Site, ncol=1)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", y="Stand LAI", alpha="Intensity")

comps_list = split(thin_swy, f=thin_swy$thin)
comps_diff = lapply(comps_list, REL_DIFF, comps_list[["Control"]], c(10:15,18,19))
comps_absdiff = lapply(comps_list, ABS_DIFF, comps_list[["Control"]], c(10:15,18,19))
stand_diff = do.call(rbind, comps_diff) %>% filter(thin!="Control")
stand_absdiff = do.call(rbind, comps_absdiff) %>% filter(thin!="Control")

ggplot(filter(stand_diff, loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+geom_hline(yintercept=0)+facet_wrap(elev~Site, ncol=1)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", y="Stand LAI (% Deviation)", alpha="Intensity")
  
ggplot(filter(stand_diff, loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=interaction(elev, Site), group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+geom_hline(yintercept=0)+facet_grid(elev~Site)+scale_fill_brewer(palette="Paired", direction=-1)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))+theme_bw()+
  theme(legend.position="bottom", axis.text=element_text(size=10), 
        strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y=expression("%"*Delta~"Stand LAI"), alpha="Intensity")+guides(fill="none")

ggplot(filter(stand_diff, loc=="Thin Both", abs(npp)<200), aes(as.factor(t_wy), npp, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+geom_hline(yintercept=0)+facet_wrap(Site~elev, ncol=1)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", y="Stand NPP (% Deviation", alpha="Intensity")

ggplot(filter(stand_absdiff, loc=="Thin Both"), aes(as.factor(t_wy), trans, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+facet_wrap(~Site, ncol=1)+geom_hline(yintercept=0)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", alpha="Intensity")

ggplot(filter(stand_diff, loc=="Thin Both"), aes(as.factor(t_wy), ET, fill=interaction(elev,Site), group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+facet_grid(elev~Site)+geom_hline(yintercept=0)+scale_fill_brewer(palette="Paired")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", alpha="Intensity")

ggplot(filter(stand_diff, loc %in% c("Thin Upslope", "Burn Upslope"), Site=="Riparian"), aes(as.factor(t_wy), ET, fill=interaction(elev,loc), group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+facet_grid(loc~elev)+geom_hline(yintercept=0)+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), 
  strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  labs(x="Year since Treatment", y=expression("%"*Delta~"ET"), alpha="Intensity")+ggtitle("Riparian Stand Response to Upslope Thin")

ggplot(filter(stand_diff, loc=="Thin Both"), aes(as.factor(t_wy), standc, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+
  geom_boxplot(outlier.size=.3)+geom_hline(yintercept=0)+facet_wrap(~Site, ncol=1)+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))+theme_bw()+theme(legend.position="bottom")+
  labs(x="Year since Treatment", y="Stand C (% Deviation)", alpha="Intensity")


# parameter uncertainty

stand_diff_p = stand_diff %>% group_by(param, Site, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% select(-thin)
CVsp_R2 = apply(subset(stand_diff_p, Site=="R", select=c(standc,lai,trans,npp)), 2, function(x){sd(x)/mean(x)})
CVsp_U2 = apply(subset(stand_diff_p, Site=="U", select=c(standc,lai,trans,npp)), 2, function(x){sd(x)/mean(x)})
CVsp2 = data.frame(Var=colnames(stand_diff_p[,8:11]), CVsp_R2, CVsp_U2)

# climate uncertainty

stand_diff_c = stand_diff %>% group_by(clim, Site, aspect, elev, int, loc, t_wy) %>% summarize_all(list(median)) %>% select(-thin)
CVsc_R2 = apply(subset(stand_diff_c, Site=="R", select=c(standc,lai,trans,npp)), 2, function(x){sd(x)/mean(x)})
CVsc_U2 = apply(subset(stand_diff_c, Site=="U", select=c(standc,lai,trans,npp)), 2, function(x){sd(x)/mean(x)})
CVs2 = cbind(CVsp2, data.frame(CVsc_R2, CVsc_U2))


CVs_comp2 = rbind(mutate(stand_diff_p, comp="Param"), mutate(stand_diff_c, comp="Clim"))
CVs_long2 = gather(CVs_comp2, key="Var", value="Val", lai, npp, standc, trans)
ggplot(CVs_long2, aes(comp, Val, fill=Site, group=interaction(comp, Site)))+geom_boxplot()+facet_wrap(~Var, scales="free_y")+
  labs(x="Uncertainty", fill="Site")


# climate comparisons

ggplot(filter(stand_diff, t_wy==2, loc=="Thin Both"), aes(pcp, trans, color=swe))+geom_point()+facet_grid(Site~elev, scales="free_y")+
  geom_hline(yintercept=0, size=.5)+labs(x="Annual P (m)", y="Annual ET (% Deviation)", color="SWE (m)")+theme_bw()+theme(strip.text=element_text(size=14))



# story metrics
thin_wy$story = with(thin_wy, case_when(patch=="011"~"Over", patch=="021"~"Over", patch=="031"~"Over", patch=="041"~"Under", patch=="051"~"Under"))
thin_wy$story_area = with(thin_wy, case_when(patch=="011"~1/7, patch=="021"~3/7, patch=="031"~3/7, patch=="041"~.5, patch=="051"~.5))
thin_sywy = thin_wy %>% group_by(Site, story, t_wy, param, clim, aspect, elev, int, loc, thin) %>% 
  summarize(standc = plantc %*% story_area, lai = lai %*% story_area, trans = trans %*% story_area, evap = evap %*% story_area, ET = ET %*% story_area, npp = npp %*% story_area, swe = swe %*% story_area, pcp = pcp %*% story_area, sddiff = sddiff %*% story_area, recharge = recharge %*% story_area, .groups="keep")
#thin_sywy$thin = if_else(thin_sywy$loc!="Control", paste(thin_sywy$loc, thin_sywy$int, sep=" "), thin_sywy$loc)

thin_story = thin_sywy %>% dplyr::select(Site:thin, ET) %>%  pivot_wider(names_from=story, values_from=ET)
thin_story$U_O = thin_story$Under / (thin_story$Under + thin_story$Over) * 100
ggplot(thin_story, aes(U_O))+geom_density()

mean(filter(thin_story, loc=="Control", Site=="Riparian")$U_O)
mean(filter(thin_story, loc!="Control", Site=="Riparian")$U_O)

compsy_list = split(thin_sywy, f=thin_sywy$thin)
compsy_diff = lapply(compsy_list, REL_DIFF, compsy_list[["Control"]], c(11:16,19,20))
compsy_absdiff = lapply(compsy_list, ABS_DIFF, compsy_list[["Control"]], c(11:16,19,20))
story_diff = do.call(rbind, compsy_diff) %>% filter(thin!="Control")
story_absdiff = do.call(rbind, compsy_absdiff) %>% filter(thin!="Control")


ggplot(filter(story_absdiff, elev=="Low-Elevation"), aes(as.factor(t_wy), ET, fill=interaction(story, Site), group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot(outlier.size=.3)+
  facet_grid(story~Site, scales="free")+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  scale_fill_brewer(palette="Paired", direction = -1)+guides(fill="none")+labs(x="Year since Treatment", y=expression(Delta~"ET"), alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

# tables

story_diff2 = story_diff %>% filter(Site=="Riparian", loc %in% c("Thin Both", "Burn Both")) %>% group_by(story) %>% 
  summarize_at(vars(lai, ET), list(mean))


# strata

ggplot(filter(thin_wy, loc %in% c("Thin Both","Control"), stratumID %in% c(38571031, 36815031, 46901031, 48536031)), aes(as.factor(t_wy), lai, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot(outlier.size=.3)+
  facet_grid(elev~strata)+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  scale_fill_brewer(palette="Paired", direction = -1)+guides(fill="none")+labs(x="Year since Treatment", y="LAI", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(thin_wy, loc=="Thin Both", stratumID %in% c(38571031, 36815031)), aes(as.factor(t_wy), height, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_wrap(~strata, ncol=1)+theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2", direction = -1)+
  labs(x="Years Since Thinning", y="Height", fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(thin_wy, loc %in% c("Thin Both","Control"), stratumID %in% c(38571041, 36815041, 46901041, 48536041)), aes(as.factor(t_wy), ET, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_grid(elev~strata)+theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Paired", direction = -1)+guides(fill="none")+
  labs(x="Year since Treatment", y="ET", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(thin_wy, loc %in% c("Thin Both","Control"), stratumID %in% c(38571051, 36815051, 46901051, 48536051)), aes(as.factor(t_wy), lai, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_grid(elev~strata)+theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Paired", direction = -1)+guides(fill="none")+
  labs(x="Year since Treatment", y="LAI", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(thin_wy, loc %in% c("Thin Both","Control"), stratumID %in% c(38571031, 36815031)), aes(as.factor(t_wy), evap, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_wrap(~strata, ncol=1)+theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2", direction = -1)+
  labs(x="Years Since Thinning", y="Evap", fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(thin_wy, loc %in% c("Thin Both","Control"), stratumID %in% c(38571051, 36815051)), aes(as.factor(t_wy), npp, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_wrap(~strata, ncol=1)+theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2", direction = -1)+
  labs(x="Years Since Thinning", y="NPP", fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


comp_long = gather(thin_wy, key="Var", value="Val", lai, npp, plantc, trans, evap, ET, height) %>% select(-c(cpool, rootdepth, sddiff, recharge, swe, pcp))
comp_long$strata = with(comp_long, case_when(stratumID==36815011~"Upslope Old-Growth", stratumID==36815021~"Ups_Med", stratumID==36815031~"Ups_Thin", stratumID==36815041~"Ups_Shrub", stratumID==36815051~"Ups_Grass",
                                           stratumID==38571011~"Riparian Old-Growth", stratumID==38571021~"Rip_Med", stratumID==38571031~"Rip_Thin", stratumID==38571041~"Rip_Shrub", stratumID==38571051~"Rip_Grass",
                                           stratumID==46901011~"Upslope Old-Growth", stratumID==46901021~"Ups_Med", stratumID==46901031~"Ups_Thin", stratumID==46901041~"Ups_Shrub", stratumID==46901051~"Ups_Grass",
                                           stratumID==48536011~"Riparian Old-Growth", stratumID==48536021~"Rip_Med", stratumID==48536031~"Rip_Thin", stratumID==48536041~"Rip_Shrub", stratumID==48536051~"Rip_Grass"))

comp_long$thin = if_else(comp_long$loc!="Control", paste(comp_long$loc, comp_long$int, sep=" "), comp_long$loc)

ggplot(filter(comp_long, Var=="lai", stratumID %in% c(38571011, 36815011)), aes(loc, Val, fill=loc))+geom_boxplot()+facet_wrap(int~strata)+
  labs(x="Treatment Type", y="Leaf Area Index (LAI)", fill="Treatment")+theme_bw()+
  theme(axis.text.x=element_text(angle=45, vjust=.5), legend.position="none")+
  scale_fill_brewer(palette="Dark2")

ggplot(filter(comp_long, Var=="lai", stratumID %in% c(38571011, 36815011, 46901011, 48536011)), aes(thin, Val, fill=thin))+geom_boxplot()+facet_wrap(elev~strata)+
  labs(x="Treatment Type", y="Leaf Area Index (LAI)", fill="Treatment")+theme_bw()+
  theme(axis.text.x=element_text(angle=45, vjust=.5), legend.position="none")+
  scale_fill_brewer(palette="Dark2")

ggplot(filter(comp_long, Var=="trans", stratumID %in% c(38571011, 36815011, 46901011, 48536011), Site=="R"), aes(as.factor(t_wy), Val, fill=loc))+geom_boxplot()+facet_wrap(~elev, ncol=1)+
  labs(y="Trans", fill="Treatment")+theme_bw()+
  theme(axis.text.x=element_text(angle=45, vjust=.5), legend.position="bottom")+
  scale_fill_brewer(palette="Dark2")

ggplot(filter(comp_long, Var=="npp", Val>-1, stratumID %in% c(38571011, 36815011, 46901011, 48536011), Site=="R"), aes(as.factor(t_wy), Val, fill=loc))+geom_boxplot()+facet_wrap(~elev, ncol=1)+
  labs(y="NPP", fill="Treatment")+theme_bw()+
  theme(axis.text.x=element_text(angle=45, vjust=.5), legend.position="bottom")+
  scale_fill_brewer(palette="Dark2")

ggplot(filter(comp_long, Var=="height", stratumID %in% c(38571011, 36815011)), aes(loc, Val, fill=loc))+geom_boxplot()+facet_wrap(int~strata)

comp_list = split(thin_wy, f=thin_wy$thin)

comp_diff = lapply(comp_list, REL_DIFF, comp_list[["Control"]], 4:14)
comp_absdiff = lapply(comp_list, ABS_DIFF, comp_list[["Control"]], 4:14)
met_diff = do.call(rbind, comp_diff) %>% filter(thin!="Control")
met_absdiff = do.call(rbind, comp_absdiff) %>% filter(thin!="Control")

ggplot(filter(met_diff, stratumID %in% c(38571011, 36815011), loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot(outlier.size=.5)+facet_wrap(~strata)+
  theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2")+
  labs(x="Years since Thinning", y="LAI (% Deviation)", fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, stratumID %in% c(38571011, 36815011, 46901011, 48536011), lai > -10, lai<50, loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)), outlier.size=.3, inherit.aes=F)+facet_grid(elev~strata, scales="free_y")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"LAI"), fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.2,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, stratumID %in% c(36815011, 48536011), lai > -5, lai<10, loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=interaction(elev, int), group=interaction(t_wy, aspect), alpha=as.factor(aspect)), outlier.size=.3, inherit.aes=F)+facet_grid(int~elev)+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Upslope Old-Growth LAI Recovery")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"LAI"), fill="Site", alpha="Aspect")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))


ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, stratumID %in% c(38571031, 36815031, 46901031, 48536031), loc=="Thin Both", lai<200), aes(as.factor(t_wy), lai, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)),outlier.size=.3)+facet_grid(elev~strata)+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+ggtitle("Thinned Tree LAI Recovery")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"LAI"), fill="Site", alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.2,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, stratumID %in% c(36815021, 48536021), lai<100, loc=="Thin Both"), aes(as.factor(t_wy), lai, fill=interaction(elev, int), group=interaction(t_wy, aspect), alpha=as.factor(aspect)), outlier.size=.3, inherit.aes=F)+facet_grid(int~elev, scales="free_y")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Upslope Medium Tree LAI Recovery")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"LAI"), fill="Site", alpha="Aspect")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))


# trans time-series
ggplot(filter(met_diff, elev=="Lo", stratumID %in% c(38571031, 36815031), trans<200, loc=="Thin Both"), aes(as.factor(t_wy), trans, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+facet_wrap(~strata, ncol=1)+
  theme_bw()+theme(legend.position="bottom")+scale_fill_brewer(palette="Dark2")+
  labs(x="Years since Thinning", y="Transpiration (% Deviation)", fill="Site")

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, loc=="Thin Both", stratumID %in% c(38571011, 36815011, 46901011, 48536011), trans > -5, trans<20), aes(as.factor(t_wy), trans, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)),outlier.shape=NA)+facet_grid(elev~strata, scales="free_y")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"Transpiration"), alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, loc=="Thin Both", stratumID %in% c(38571011, 36815011, 46901011, 48536011), trans > -5, trans<15), aes(aspect, trans, fill=interaction(elev, strata), group=interaction(aspect, int), alpha=as.factor(int)),outlier.shape=NA)+facet_grid(elev~strata)+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+
  labs(x="Aspect", y=expression("%"*Delta~"Transpiration"), alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, loc=="Thin Both", stratumID %in% c(38571011, 36815011, 46901011, 48536011)), aes(as.factor(t_wy), trans, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)),outlier.size=.3)+facet_grid(elev~strata, scales="free_y")+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"Transpiration"), alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, loc=="Thin Both", stratumID %in% c(38571011, 36815011), trans > -25, trans<60), aes(as.factor(t_wy), trans, fill=interaction(int, strata), group=interaction(t_wy, aspect), alpha=as.factor(aspect)), outlier.size=.3)+facet_grid(int~strata)+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Year since Treatment", y=expression("%"*Delta~"Transpiration"), alpha="Aspect")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_diff, loc=="Thin Both", trans<100), aes(patch, trans, fill=interaction(elev, patch), group=interaction(int, patch), alpha=as.factor(int)), outlier.size=.3)+facet_grid(elev~Site)+
  theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Dark2")+guides(fill="none")+
  labs(x="Canopy Layer", y=expression("%"*Delta~"Transpiration"), alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))


ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_absdiff, loc %in% c("Thin Both", "Burn Both"), stratumID %in% c(38571011, 36815011, 46901011, 48536011)), aes(as.factor(t_wy), ET, fill=interaction(elev, strata), group=interaction(t_wy, int), alpha=as.factor(int)), outlier.size=.3)+
  facet_grid(elev~strata, scales="free_y")+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+scale_fill_brewer(palette="Paired", direction=-1)+guides(fill="none")+
  labs(x="Year since Treatment", y=expression(Delta~"ET"), alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_absdiff, stratumID %in% c(38571041, 48536041), loc %in% c("Thin Both", "Burn Both")), aes(as.factor(t_wy), ET, fill=interaction(elev, loc), group=interaction(t_wy, aspect), alpha=as.factor(aspect)), outlier.size=.3, inherit.aes=F)+
  facet_grid(loc~elev)+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Riparian Shrub ET Recovery")+
  labs(x="Year since Treatment", y=expression(Delta~"ET (mm/year)"), alpha="Aspect")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_absdiff, stratumID %in% c(36815041, 46901041), loc %in% c("Thin Both", "Burn Both")), aes(as.factor(t_wy), ET, fill=interaction(elev, loc), group=interaction(t_wy, aspect), alpha=as.factor(aspect)), outlier.size=.3, inherit.aes=F)+
  facet_grid(loc~elev)+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Upslope Shrub ET Recovery")+
  labs(x="Year since Treatment", y=expression(Delta~"ET (mm/year)"), alpha="Aspect")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_absdiff, stratumID %in% c(36815051, 46901051), loc %in% c("Thin Both", "Burn Both")), aes(as.factor(t_wy), ET, fill=interaction(elev, loc), group=interaction(t_wy, int), alpha=as.factor(int)), outlier.size=.3, inherit.aes=F)+
  facet_grid(loc~elev)+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Upslope Grass ET Recovery")+
  labs(x="Year since Treatment", y=expression(Delta~"ET (mm/year)"), alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))

ggplot()+geom_hline(yintercept=0)+
  geom_boxplot(data=filter(met_absdiff, stratumID %in% c(38571051, 48536051), loc %in% c("Thin Both", "Burn Both")), aes(as.factor(t_wy), ET, fill=interaction(elev, loc), group=interaction(t_wy, int), alpha=as.factor(int)), outlier.size=.3, inherit.aes=F)+
  facet_grid(loc~elev)+theme_bw()+theme(legend.position="bottom", axis.text=element_text(size=10), strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  scale_fill_brewer(palette="Dark2")+guides(fill="none")+ggtitle("Riparian Grass ET Recovery")+
  labs(x="Year since Treatment", y=expression(Delta~"ET (mm/year)"), alpha="Intensity")+scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")), range=c(.3,.9))


ggplot(filter(met_diff, loc %in% c("Thin Both", "Burn Both"), stratumID %in% c(38571041, 36815041)), aes(as.factor(t_wy), ET, fill=Site, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_grid(loc~strata)+theme_bw()+theme(legend.position="bottom")+
  labs(x="Years since Thinning", y="ET (% Deviation)", fill="Location", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))

ggplot(filter(met_diff, loc %in% c("Thin Upslope", "Burn Upslope"), stratumID %in% c(38571031, 48536031)), aes(as.factor(t_wy), ET, fill=strata, group=interaction(t_wy, int), alpha=as.factor(int)))+geom_boxplot()+
  facet_grid(loc~elev)+theme_bw()+theme(legend.position="bottom")+
  labs(x="Years since Thinning", y="ET (% Deviation)", fill="Location", alpha="Intensity")+
  scale_alpha_discrete(guide = guide_legend(override.aes = list(fill = "black")))


# stand scale water balance

met_wy = met_absdiff %>% select(-c(plantc, swe, pcp, recharge, cpool, height, rootdepth)) %>% 
  group_by(t_wy, Site, patch, elev, int, loc, thin) %>% summarize_all(list(mean)) %>% 
  select(-c(scen, param, clim, aspect, type, strata, stratumID))

met_wy2 = met_diff %>% select(-c(plantc, swe, pcp, recharge, cpool, height, rootdepth)) %>% 
  group_by(t_wy, Site, patch, strata, stratumID, elev, int, loc, thin) %>% summarize_all(list(mean)) %>% 
  select(-c(scen, param, clim, aspect, type))

ggplot(filter(met_wy, elev=="Low-Elevation", thin %in% c("Burn Both 30","Thin Both 30"), patch != "031"), aes(t_wy, ET, fill=patch))+
  geom_col()+facet_grid(factor(thin, levels=c("Thin Both 30","Burn Both 30"))~Site, scales="free")+theme_bw()+
  theme(legend.position="bottom", axis.text=element_text(size=10), 
        strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"), plot.title=element_text(hjust=.5))+
  labs(x="Year since Treatment", y=expression(Delta~"ET (mm/yr)"), fill="Canopy Layer")+
  scale_fill_discrete(labels=c("Old-Growth","Medium Trees", "Shrub", "Grass"))+
  ggtitle("Comparison of Mean Change in ET from Treatment by Canopy Layer")

ggplot(filter(met_wy, elev=="Low-Elevation", thin %in% c("Burn Both 30","Thin Both 30")), aes(t_wy, lai, color=patch))+
  geom_line(linewidth=1)+facet_grid(Site~thin)

ggplot(filter(met_wy2, elev=="Low-Elevation", thin %in% c("Burn Both 30","Thin Both 30")), aes(t_wy, npp, color=patch))+
  geom_line()+facet_grid(factor(thin, levels=c("Thin Both 30","Burn Both 30"))~Site)+theme_bw()+
  theme(legend.position="bottom", axis.text=element_text(size=10), 
        strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))+
  labs(x="Year since Treatment", y=expression("%"*Delta~"NPP"), color="Canopy Layer")+
  scale_color_discrete(labels=c("Old-Growth Trees","Medium Trees", "Thinned Trees", "Shrub", "Grass"))


met_wb = met_wy %>% select(-c(trans, evap, lai, npp)) %>% rename(S=sddiff) %>% 
  pivot_wider(names_from=patch, values_from=c(ET, S))
met_wb2 = gather(met_wb, key="Flux", value="Value", ET_011:S_051)
met_wb2$Flux = factor(met_wb2$Flux, levels=c("S_011","ET_011","S_021","ET_021","S_031","ET_031","S_041","ET_041","S_051","ET_051"),
                      labels=c("Old-Growth Storage","Old-Growth ET","Medium Tree Storage","Medium ET",
                               "Thinned Tree Storage","Thinned Tree ET","Shrub Storage","Shrub ET","Grass Storage","Grass ET"))

ggplot(filter(met_wb2, elev=="Low-Elevation", thin=="Burn Both 30"), aes(t_wy, Value, fill=Flux))+
  geom_col()+facet_wrap(~Site, ncol=1)+theme_bw()+labs(x="Year Since Treatment", y="Change in Water Balance Flux")+
  scale_fill_brewer(palette="Paired")+theme(legend.position="bottom")

ggplot(filter(met_wb2, elev=="Low-Elevation", thin %in% c("Burn Both 30","Thin Both 30")), aes(t_wy, Value, color=Flux))+
  geom_line(linewidth=1)+facet_grid(factor(thin, levels=c("Thin Both 30", "Burn Both 30"))~Site)+theme_bw()+labs(x="Year Since Treatment", y=expression(Delta~"Water Balance Flux (mm/yr)"))+
  scale_color_brewer(palette="Paired")+theme(legend.position="bottom", axis.text=element_text(size=10), 
                                             strip.text.x=element_text(size=10, face="bold"),strip.text.y=element_text(size=10, face="bold"))